{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3652aca5-46be-454d-9b5e-5b795379c534",
   "metadata": {},
   "source": [
    "# Modeling an HHL Algorithm to Solve a Set of Linear Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "Guidance for the workshop:\n",
    "**The `# TODO` is there for you to do yourself.**\n",
    "**The `# Solution start` and `# Solution end` are only for helping you. Try doing it yourself...**"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "id": "422c4e26-b10c-4d77-b7ea-5d06d684a752",
   "metadata": {},
   "source": [
    "Solving linear equations appears in many research, engineering, and design fields. For example, many physical and financial models, from fluid dynamics to portfolio optimization, are described by partial differential equations, which are typically treated by numerical schemes, most of which are eventually transformed to a set of linear equations.\n",
    "\n",
    "The HHL algorithm [[1](#HHL)] is a quantum algorithm for solving a set of linear equations. It is one of the fundamental quantum algorithms that is expected to give a speedup over its classical counterpart.\n",
    "\n",
    "\n",
    "A set of linear equations of size $N$ is represented by an $N\\times N$ matrix and a vector $b$ of size $N$, $A\\vec{x} = \\vec{b}$, where the solution to the problem is designated by the solution variable $\\vec{x}$.\n",
    "\n",
    "For simplicity, the demo below treats a usecase where $\\vec{b}$ is a normalized vector $|\\vec{b}|=1$, and $A$ is an Hermitian matrix of size $2^n\\times 2^n$, whose eigenvalues are in the interval $(0,1)$. Generalizations to other usecases are discussed at the end of this demo."
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 1. Defining a Specific Problem"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "id": "fb8942f5-23a0-4370-b83b-d6c41f6af033",
   "metadata": {},
   "source": [
    "Start by defining the specific problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# !pip install -U classiq"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as scipy\n",
    "\n",
    "A = np.array(\n",
    "    [\n",
    "        [0.28, -0.01, 0.02, -0.1],\n",
    "        [-0.01, 0.5, -0.22, -0.07],\n",
    "        [0.02, -0.22, 0.43, -0.05],\n",
    "        [-0.1, -0.07, -0.05, 0.42],\n",
    "    ]\n",
    ")\n",
    "\n",
    "b = np.array([1, 2, 4, 3])\n",
    "b = b / np.linalg.norm(b)\n",
    "\n",
    "print(\"A =\", A, \"\\n\")\n",
    "print(\"b =\", b, \"\\n\")\n",
    "\n",
    "# verifying that the matrix is symmetric and has eigenvalues in [0,1)\n",
    "if not np.allclose(A, A.T, rtol=1e-6, atol=1e-6):\n",
    "    raise Exception(\"The matrix is not symmetric\")\n",
    "w, v = np.linalg.eig(A)\n",
    "for lam in w:\n",
    "    if lam < 0 or lam > 1:\n",
    "        raise Exception(\"Eigenvalues are not in (0,1)\")\n",
    "\n",
    "sol_classical = np.linalg.solve(A, b)\n",
    "print(\"Classical solution: x = \", sol_classical)\n",
    "num_qubits = int(np.log2(len(b)))"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 2. Building Simple HHL with Classiq"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2.1 Define The Model\n",
    "This tutorial gives instructions on building an HHL algorithm and presents the theory of the algorithm. The algorithm consists of 4 steps:\n",
    "\n",
    "1) State preparation of the RHS vector $\\vec{b}$.\n",
    "\n",
    "2) QPE for the unitary matrix $e^{2\\pi iA}$, which encodes eigenvalues on a quantum register of size $m$.\n",
    "\n",
    "3) An inversion algorithm that loads amplitudes according to the inverse of the eigenvalue registers.\n",
    "\n",
    "4) An inverse QPE with the parameters in (2)."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 2.1.1 State Preparation for the Vector $\\vec{b}$\n",
    "\n",
    "The first stage of the HHL algorithm is to load the normalized RHS vector $\\vec{b}$ into a quantum register:\n",
    "\n",
    "$$\n",
    "|0\\rangle_n \\xrightarrow[{\\rm SP}]{} \\sum^{2^n-1}_{i=0}b_i|i\\rangle_n\n",
    "$$\n",
    "\n",
    "where $|i\\rangle$ are states in the computational basis.\n",
    "\n",
    "Comments:\n",
    "\n",
    "* The relevant built-in function is the `prepare_amplitudes` one, which gets $2^n$ values of $\\vec{b}$, as well as an upper bound for its functional error through the `bound` parameter."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import CArray, CReal, Output, QArray, QBit, prepare_amplitudes, qfunc\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def load_b(b: CArray[CReal], res: Output[QArray[QBit]]) -> None:\n",
    "    # TODO prepare the state |b> in the \"res\" register - the amplitude of res states correspond to the values of the vector b\n",
    "    # Solution start\n",
    "    prepare_amplitudes(b, 0.0, res)\n",
    "    # Solution end"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's see the loading of b in a state vector simulator\n",
    "\n",
    "Update the qmod to have `aer_simulator_statevector` as backend, with one shot\n",
    "\n",
    "Refer to [Execution Preferences documentation](https://docs.classiq.io/latest/reference-manual/platform/executor/#execution-preferences) and to [Classiq backends documentation](https://docs.classiq.io/latest/reference-manual/platform/executor/cloud-providers/)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def main(res: Output[QArray[QBit]]):\n",
    "    load_b(b.tolist(), res)\n",
    "\n",
    "\n",
    "from classiq import create_model, execute, show, synthesize\n",
    "from classiq.execution import ClassiqBackendPreferences, ExecutionPreferences\n",
    "from classiq.synthesis import set_execution_preferences\n",
    "\n",
    "qmod_b_load = create_model(main)\n",
    "# TODO update the qmod to have aer_simulator_statevector as backend, with one shot\n",
    "\n",
    "# Solution start\n",
    "backend_preferences = ClassiqBackendPreferences(\n",
    "    backend_name=\"aer_simulator_statevector\"\n",
    ")\n",
    "qmod_b_load = set_execution_preferences(\n",
    "    qmod_b_load,\n",
    "    execution_preferences=ExecutionPreferences(\n",
    "        num_shots=1, backend_preferences=backend_preferences\n",
    "    ),\n",
    ")\n",
    "# Solution end\n",
    "qprog_b_load = synthesize(qmod_b_load)\n",
    "show(qprog_b_load)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Take a look at the resulted circuit.\n",
    "Now let's execute and see if $b$ was built correctly"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "job = execute(qprog_b_load)\n",
    "job.open_in_ide()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Check if you see a match between the original $b$ to the resulted state vector"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "print(\"The original b is: \", b)\n",
    "results = job.result()\n",
    "res_hhl = results[0].value\n",
    "print(\"The resulted state vector :\", res_hhl.state_vector)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 2.1.2 Quantum Phase Estimation (QPE) for the Hamiltonian Evolution $U=e^{2\\pi i A }$\n",
    "\n",
    "The QPE function block, which is at the heart of the HHL algorithm, operates as follows: Unitary matrices have eigenvalues of norm 1, and thus are of the form $e^{2\\pi i \\lambda}$, with $0\\leq\\lambda<1$. For a quantum state $|\\psi\\rangle_n$, prepared in an eigenvalue of some unitary matrix $U$ of size $2^n\\times 2^n$, the QPE algorithm encodes the corresponding  eigenvalue into a quantum register:\n",
    "\n",
    "$$\n",
    "|0\\rangle_m|\\psi\\rangle_n \\xrightarrow[{\\rm QPE}(U)]{} |\\lambda\\rangle_m|\\psi\\rangle_n,\n",
    "$$\n",
    "\n",
    "where $m$ is the precision of the binary representation of $\\lambda$, $\\lambda=\\frac{1}{2^m}\\sum^{2^m-1}_{k=0}\\lambda^{(k)}2^k$ with $\\lambda^{(k)}$ being the state of the $k$-th qubit.\n",
    "\n",
    "In the HHL algorithm a QPE for the unitary $U=e^{2\\pi i A }$ is applied. The mathematics: First, note that the eigenvectors of $U$ are the ones of the matrix $A$, and that the corresponding eigenvalues $\\lambda$ defined for $U=e^{2\\pi i A }$ are the eigenvalues of $A$. Second, represent the prepared state in the basis given by the eigenvalues of $A$. This is merely a mathematical transformation; with no algorithmic considerations here. If the eigenbasis of $A$ is given by the set $\\{|\\psi_j\\rangle_n \\}^{2^n-1}_{j=0}$, then\n",
    "\n",
    "$$\n",
    "\\sum^{2^n-1}_{i=0}b_i|i\\rangle_n = \\sum^{2^n-1}_{j=0}\\beta_j|\\psi_j\\rangle_n.\n",
    "$$\n",
    "\n",
    "Applying the QPE stage gives\n",
    "\n",
    "$$\n",
    "|0\\rangle_m \\sum^{2^n-1}_{j=0}\\beta_j|\\psi_j\\rangle_n \\xrightarrow[{\\rm QPE}]{}  \\sum^{2^n-1}_{j=0}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n.\n",
    "$$\n",
    "\n",
    "Comments:\n",
    "\n",
    "* Use the built-in `qpe` function."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 2.1.3 Eigenvalue Inversion\n",
    "The next step in the HHL algorithm is to pass the inverse of the eigenvalue registers into their amplitudes, using the Amplitude Loading (AL) construct.\n",
    "Given a function $f:[0,1)\\rightarrow [-1,1]$, it implements $|0\\rangle|\\lambda\\rangle_m \\xrightarrow[{\\rm AL}(f)]{} f(\\lambda)|1\\rangle|\\lambda\\rangle_m+\\sqrt{1-f^2(\\lambda)}|0\\rangle|\\lambda\\rangle_m$. For the HHL algorithm, apply an AL with $f=C/x$ where $C$ is a lower bound for the minimal eigenvalue of $A$.\n",
    "Applying this AL gives\n",
    "\n",
    "$$\n",
    "\\sum^{2^n-1}_{j=0}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n \\xrightarrow[{\\rm AL}(C/x)]{}\n",
    "|0\\rangle\\left(\\sum^{2^n-1}_{j=0}\\sqrt{1-\\frac{C^2}{\\lambda^2_j}}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n\\right)+\n",
    "|1\\rangle\\left(\\sum^{2^n-1}_{j=0}\\frac{C}{\\lambda_j}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n\\right),\n",
    "$$\n",
    "\n",
    "where $C$ is a normalization coefficient.\n",
    "\n",
    "The normalization coefficient $C$, which guarantees that the amplitudes are normalized, can be taken as the lower possible eigenvalue that can be resolved with the QPE:\n",
    "\n",
    "$$\n",
    "C=1/2^{\\rm precision}.\n",
    "$$\n",
    "\n",
    "The built-in construct to define an amplitude loading is designated by `*=`, see [Expression Assignment Operations](https://docs.classiq.io/latest/reference-manual/platform/qmod/language-reference/statements/numeric-assignment/#example-3-in-place-assignment-of-a-logical-expression)."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import QNum, allocate\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def simple_eig_inv(phase: QNum, indicator: Output[QBit]):\n",
    "    # TODO allocate 1 qubit for indicator\n",
    "    # TODO load its |1> state amplitude to be C/phase using the *= operator\n",
    "\n",
    "    # Solution start\n",
    "    allocate(1, indicator)\n",
    "    indicator *= (1 / 2**phase.size) / phase\n",
    "    # Solution end"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 2.1.4 Inverse QPE\n",
    "\n",
    "As the final step in the HHL model, clean the QPE register by applying an inverse-QPE. (Note that it is not guaranteed that this register is completely cleaned; namely, that all the qubits in the QPE register return to zero after the inverse-QPE. Generically they are all zero with very high probability).\n",
    "\n",
    "In this model we will simply call the QPE function with the same parameters in stage 2. This is how the quantum state looks now\n",
    "\n",
    "$$\n",
    " |0\\rangle\\left(\\sum^{2^n-1}_{j=0}\\sqrt{1-\\frac{C^2}{\\lambda^2_j}}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n\\right)+\n",
    "|1\\rangle\\left(\\sum^{2^n-1}_{j=0}\\frac{C}{\\lambda_j}\\beta_j |\\lambda_j\\rangle_m |\\psi_j\\rangle_n\\right)\n",
    "\\xrightarrow[{\\rm inv-QPE}(U)]{}\n",
    " |0\\rangle_m|0\\rangle\\left(\\sum^{2^n-1}_{j=0}\\sqrt{1-\\frac{C^2}{\\lambda^2_j}}\\beta_j  |\\psi_j\\rangle_n\\right)+\n",
    "|0\\rangle_m|1\\rangle\\left(\\sum^{2^n-1}_{j=0}\\frac{C}{\\lambda_j}\\beta_j  |\\psi_j\\rangle_n\\right)\n",
    "$$\n",
    "\n",
    "The state entangled with $|1\\rangle$ stores the solution to our problem (up to some normalization problem)\n",
    "\n",
    "$$\n",
    "\\sum^{2^n-1}_{j=0} \\frac{C}{\\lambda_j}\\beta_j  \\vec{\\psi_j} = C\\vec{x}.\n",
    "$$\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 2.1.5 Putting it all together\n",
    "\n",
    "Let's remind that the entire HHL algorithm is composed of:\n",
    "\n",
    "1) State preparation of the RHS vector $\\vec{b}$.\n",
    "\n",
    "2) QPE for the unitary matrix $e^{2\\pi iA}$, which encodes eigenvalues on a quantum register of size $m$.\n",
    "\n",
    "3) An inversion algorithm that loads amplitudes according to the inverse of the eigenvalue registers.\n",
    "\n",
    "4) An inverse QPE with the parameters in (2).\n",
    "\n",
    "And put all together in `my_hhl` function\n",
    "\n",
    "You can apply QPE† * EigenValInv * QPE using the [within_apply operator](https://docs.classiq.io/latest/reference-manual/platform/qmod/language-reference/statements/within-apply/)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import CInt, QCallable, allocate_num, qpe, unitary, within_apply\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def my_hhl(\n",
    "    fraction_digits: CInt,\n",
    "    b: CArray[CReal],\n",
    "    unitary_with_matrix: QCallable[QArray[QBit]],\n",
    "    res: Output[QArray[QBit]],\n",
    "    phase: Output[QNum],\n",
    "    indicator: Output[QBit],\n",
    ") -> None:\n",
    "    # TODO Call load_b you created, to load \"b\" vector into register \"res\"\n",
    "    # Solution start\n",
    "    load_b(b, res)\n",
    "    # Solution end\n",
    "\n",
    "    # TODO allocate a qnum register for \"phase\". This qnum should be in the range [0,1) with fraction_digits precision\n",
    "    # Solution start\n",
    "    allocate_num(fraction_digits, False, fraction_digits, phase)\n",
    "    # Solution end\n",
    "\n",
    "    # TODO refer to applying (QPE†)*(EigenValInv)*(QPE) : we want to apply \"simple_eig_inv\" within \"qpe\"\n",
    "    # Solution start\n",
    "    within_apply(\n",
    "        lambda: qpe(unitary=lambda: unitary_with_matrix(res), phase=phase),\n",
    "        lambda: simple_eig_inv(phase=phase, indicator=indicator),\n",
    "    )\n",
    "    # Solution end"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "The first entry point of any model would be the `main` function.\n",
    "Since you already have done all the job in `my_hhl` function, all we have to do now is to call it with the relevant inputs\n",
    "\n",
    "Call `my_hhl` with `QPE_SIZE` digits resolution, on the normalized `b`, where the unitary is based on `unitary_mat`."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "QPE_SIZE = 4\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def main(res: Output[QArray[QBit]], phase: Output[QNum], indicator: Output[QBit]):\n",
    "    b_normalized = b.tolist()\n",
    "    unitary_mat = scipy.linalg.expm(1j * 2 * np.pi * A).tolist()\n",
    "\n",
    "    # TODO call my_hhl with QPE_SIZE digits resolution, on the normalized b, where the unitary is based on \"unitary_mat\"\n",
    "\n",
    "    # Solution start\n",
    "    my_hhl(\n",
    "        fraction_digits=QPE_SIZE,\n",
    "        b=b_normalized,\n",
    "        unitary_with_matrix=lambda target: unitary(elements=unitary_mat, target=target),\n",
    "        res=res,\n",
    "        phase=phase,\n",
    "        indicator=indicator,\n",
    "    )\n",
    "    # Solution end"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2.2 Add Execution Preferences\n",
    "Once we have a model `qmod_hhl` (by creating it from `main`), we would like to add execution details to prepare for the program execution stage."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import show, synthesize\n",
    "from classiq.execution import ClassiqBackendPreferences, ExecutionPreferences\n",
    "from classiq.synthesis import set_execution_preferences\n",
    "\n",
    "backend_preferences = ClassiqBackendPreferences(\n",
    "    backend_name=\"aer_simulator_statevector\"\n",
    ")\n",
    "\n",
    "qmod_hhl = create_model(\n",
    "    entry_point=main,\n",
    "    execution_preferences=ExecutionPreferences(\n",
    "        num_shots=1, backend_preferences=backend_preferences\n",
    "    ),\n",
    ")"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2.2 Synthesize - from qmod to qprog\n",
    "Once we have a high level model, we would like to compile and get the actual quantum program.\n",
    "This is done using the `synthesize` command."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "qprog_hhl = synthesize(qmod_hhl)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Viewing in IDE"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "show(qprog_hhl)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Details about `qprog` as depth for example, can be seen both in IDE and in Python SDK"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import QuantumProgram\n",
    "\n",
    "circuit_hhl = QuantumProgram.from_qprog(qprog_hhl)\n",
    "print(\"depth = \", circuit_hhl.transpiled_circuit.depth)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "print(\"depth = \", circuit_hhl.transpiled_circuit.get_circuit_metrics())"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2.3 Execution\n",
    "Here we execute the circuit on state vector simulator (the backend for execution is defined before the synthesis stage).\n",
    "We can show the results in the IDE and save the state-vector into a variable."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from classiq import execute\n",
    "\n",
    "job = execute(qprog_hhl)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "job.open_in_ide()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "results = job.result()\n",
    "res_hhl = results[0].value"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2.4 Post-process\n",
    "We would like to look at the answers that are encoded in the amplitudes of the terms that their indicator qubit value=1"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "target_pos = res_hhl.physical_qubits_map[\"indicator\"][0]  # position of control qubit\n",
    "sol_pos = list(res_hhl.physical_qubits_map[\"res\"])  # position of solution\n",
    "phase_pos = list(\n",
    "    res_hhl.physical_qubits_map[\"phase\"]\n",
    ")  # position of the “phase” register, and flips for endianness as we will use the indices to read directly from the string"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Define a run over all the relevant strings holding the solution. The solution vector will be inserted into the variable `qsol`. Factor out $C=1/2^m$."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "qsol = [\n",
    "    np.round(parsed_state.amplitude / (1 / 2**QPE_SIZE), 5)\n",
    "    for solution in range(2**num_qubits)\n",
    "    for parsed_state in res_hhl.parsed_state_vector\n",
    "    if parsed_state[\"indicator\"] == 1.0\n",
    "    and parsed_state[\"res\"] == solution\n",
    "    and parsed_state[\"phase\"]\n",
    "    == 0.0  # this takes the entries where the “phase” register is at state zero\n",
    "]\n",
    "print(\"Quantum Solution: \", np.abs(qsol) / np.linalg.norm(qsol))\n",
    "print(\"Classical Solution: \", sol_classical / np.linalg.norm(sol_classical))"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "fidelity = (\n",
    "    np.abs(\n",
    "        np.dot(\n",
    "            sol_classical / np.linalg.norm(sol_classical),\n",
    "            qsol / np.linalg.norm(qsol),\n",
    "        )\n",
    "    )\n",
    "    ** 2\n",
    ")\n",
    "print(\"Solution Fidelity:\", fidelity)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# 3. Comparing classical and quantum solutions.\n",
    "\n",
    "Note that the HHL algorithm returns a statevector result up to some global phase (coming from transpilation or from the quantum functions themselves). Therefore, to compare with the classical solution, correct for this global phase.\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "sol_classical = np.linalg.solve(A, b)\n",
    "global_phase = np.angle(qsol)\n",
    "qsol_corrected = np.real(qsol / np.exp(1j * global_phase))\n",
    "print(\"classical:  \", sol_classical)\n",
    "print(\"HHL:        \", qsol_corrected)\n",
    "print(\n",
    "    \"relative distance:  \",\n",
    "    round(\n",
    "        np.linalg.norm(sol_classical - qsol_corrected)\n",
    "        / np.linalg.norm(sol_classical)\n",
    "        * 100,\n",
    "        1,\n",
    "    ),\n",
    "    \"%\",\n",
    ")"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(sol_classical, \"bo\", label=\"classical\")\n",
    "plt.plot(qsol_corrected, \"ro\", label=\"HHL\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"$i$\")\n",
    "plt.ylabel(\"$x_i$\")\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "id": "c5d5ea31-3ab2-4eb5-8522-9462e6ccfa40",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 4. Generalizations\n",
    "\n",
    "The usecase treated above is a canonical one, assuming the following properties:\n",
    "\n",
    "1) The RHS vector $\\vec{b}$ is normalized.\n",
    "\n",
    "2) The matrix $A$ is an Hermitian one.\n",
    "\n",
    "3) The matrix $A$ is of size $2^n\\times 2^n $.\n",
    "\n",
    "4) The eigenvalues of $A$ are in the range $(0,1)$.\n",
    "\n",
    "However, any general problem that does not follow these conditions can be resolved as follows:\n",
    "\n",
    "1) As preprocessing, normalize $\\vec{b}$ and then return the normalization factor as a post-processing\n",
    "\n",
    "2) Symmetrize the problem as follows:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "0 & A^T \\\\\n",
    "A & 0\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\vec{b}  \\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "0  \\\\\n",
    "\\vec{x}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "This increases the number of qubits by 1.\n",
    "\n",
    "3) Complete the matrix dimension to the closest $2^n$ with an identity matrix. The vector $\\vec{b}$ will be completed with zeros.\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "A & 0 \\\\\n",
    "0 & I\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\vec{b}  \\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "\\vec{x}  \\\\\n",
    "0\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "4) If the eigenvalues of $A$ are in the range $[-w_{\\min},w_{\\max}]$ you can employ transformations to the exponentiated matrix that enters into the Hamiltonian simulation, and then undo them for extracting the results:\n",
    "\n",
    "$$\n",
    "\\tilde{A}=(A+w_{\\min}I)\\left(1-\\frac{1}{2^{m}}\\right)\\frac{1}{w_{\\min}+w_{\\max}}.\n",
    "$$\n",
    "\n",
    "The eigenvalues of this matrix lie in the interval $[0,1)$, and are related to the eigenvalues of the original matrix via\n",
    "\n",
    "$$\n",
    "\\lambda = (w_{\\min}+w_{\\max})\\tilde{\\lambda}\\left[1/\\left(1-\\frac{1}{2^{n_{m}}}\\right)\\right]-w_{\\min},\n",
    "$$\n",
    "\n",
    "with $\\tilde{\\lambda}$ being an eigenvalue of $\\tilde{A}$ resulting from the QPE algorithm. This relation between eigenvalues is then used for the expression inserted into the eigenvalue inversion, via the `AmplitudeLoading` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfa04f60-9f81-49d9-860e-179169bd5cdf",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a name='HHL'>[1]</a>: [Harrow, A. W., Hassidim, A., & Lloyd, S., Quantum Algorithm for Linear Systems of Equations. Physical Review Letters 103, 150502 (2009)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.103.150502).\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
